#ifndef _WFIP1_H_
#define _WFIP1_H_

#include "./TwoPhaseEquation.hpp"

namespace Tuna {

  template <typename T, int Dim> class WFIP1;

  template <typename T>
  struct Typeinfo <WFIP1<T, 1> > {
    typedef T prec_t;
    enum { Dim = 1 };
  };

  template <typename T>
  struct Typeinfo <WFIP1<T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };

  template <typename T>
  struct Typeinfo <WFIP1<T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };

  
  template<typename Tprec, int Dim>
  class WFIP1 : public TwoPhaseEquation<WFIP1<Tprec, Dim> >
  {

    typedef TwoPhaseEquation<WFIP1<Tprec, Dim> > TP_WFIP1;
    
    using GeneralEquation< TP_WFIP1 >::phi;
    using GeneralEquation< TP_WFIP1 >::aE;
    using GeneralEquation< TP_WFIP1 >::aW;
    using GeneralEquation< TP_WFIP1 >::aN;
    using GeneralEquation< TP_WFIP1 >::aS;
    using GeneralEquation< TP_WFIP1 >::aB;
    using GeneralEquation< TP_WFIP1 >::aP;
    using GeneralEquation< TP_WFIP1 >::sp;
    using GeneralEquation< TP_WFIP1 >::dx;
    using GeneralEquation< TP_WFIP1 >::dy;
    using GeneralEquation< TP_WFIP1 >::dz;
    using GeneralEquation< TP_WFIP1 >::bi;
    using GeneralEquation< TP_WFIP1 >::ei;
    using GeneralEquation< TP_WFIP1 >::bj;
    using GeneralEquation< TP_WFIP1 >::ej;
    using GeneralEquation< TP_WFIP1 >::bk;
    using GeneralEquation< TP_WFIP1 >::ek;
    using GeneralEquation< TP_WFIP1 >::applyBoundaryConditions1D;
    using GeneralEquation< TP_WFIP1 >::applyBoundaryConditions2D;
    using GeneralEquation< TP_WFIP1 >::applyBoundaryConditions3D;

    using TP_WFIP1::S;
    //  using TP_WFIP1::p;
    using TP_WFIP1::phi_0;
    using TP_WFIP1::Srw;
    using TP_WFIP1::Sro;
    using TP_WFIP1::mu_w;
    using TP_WFIP1::mu_o;
    using TP_WFIP1::k;
    using TP_WFIP1::theta;
    using TP_WFIP1::p_threshold;

  public:
    typedef Tprec prec_t;
    typedef typename TunaArray<prec_t, Dim >::huge ScalarField;
    
    WFIP1() : TwoPhaseEquation<WFIP1<prec_t, Dim > >() { }    
    ~WFIP1() { };
    
    inline bool calcCoefficients1D(); 
    inline bool calcCoefficients2D() { };
    inline bool calcCoefficients3D() { };
    inline void printInfo() { std::cout << " WFIP1 "; }
  };

  
/*
 *  Lineal for realtive permeability , Upwind for Sw
 */
template<typename Tprec, int Dim>
inline bool WFIP1<Tprec, Dim>::calcCoefficients1D () 
{
    static prec_t Sw_e, Sw_w;
    static prec_t k_dx_mu_o = k / (dx * mu_o);
    static prec_t k_dx_mu_w = k / (dx * mu_w);
    static prec_t SrwSro1 = 1 - Srw - Sro;
    static prec_t SrwSro2 = SrwSro1 * SrwSro1;
    static prec_t SrwSro4 = SrwSro2 * SrwSro2;
    static prec_t SrwSro7 = SrwSro2 * SrwSro2 * SrwSro2 * SrwSro1;
    static prec_t COEF1 = k_dx_mu_w / SrwSro4;
    static prec_t COEF2 = k_dx_mu_o / SrwSro2;
    static prec_t COEF3 = k_dx_mu_o / ( SrwSro2 * SrwSro4 );
    static prec_t Sw_eSrw, Sw_eSrw4, Sw_wSrw, Sw_wSrw4;
    static prec_t SroSw_e2, SroSw_w2;
    static prec_t COEF4 = -k * p_threshold / 
      (mu_w * theta * dx * sqrt(SrwSro7));
    static prec_t aPs, aEs, aWs;

    aE = 0.0; aW = 0.0; aP = 0.0; sp = 0.0;

    for (int i =  bi; i <= ei; ++i) {  

      // Upwind
      if ( phi_0(i+1) >= phi_0(i) ) Sw_e = S(i+1);
      else                          Sw_e = S(i);
      if ( phi_0(i-1) >= phi_0(i) ) Sw_w = S(i-1);
      else                          Sw_w = S(i);
      
      Sw_eSrw  = Sw_e - Srw;
      Sw_eSrw4 = Sw_eSrw * Sw_eSrw * Sw_eSrw * Sw_eSrw;

      Sw_wSrw  = Sw_w - Srw;
      Sw_wSrw4 = Sw_wSrw * Sw_wSrw * Sw_wSrw * Sw_wSrw;

      SroSw_e2 = (1 - Sro - Sw_e) * (1 - Sro - Sw_e);
      SroSw_w2 = (1 - Sro - Sw_w) * (1 - Sro - Sw_w);

      aE (i) = COEF1 * Sw_eSrw4 + COEF2 * SroSw_e2 
	     - COEF3 * SroSw_e2 * Sw_eSrw4;
      aW (i) = COEF1 * Sw_wSrw4 + COEF2 * SroSw_w2 
	     - COEF3 * SroSw_w2 * Sw_wSrw4;
      aP (i) = aE (i) + aW (i);
      
      aEs = COEF4 * sqrt( Sw_eSrw4 * Sw_eSrw );
      aWs = COEF4 * sqrt( Sw_wSrw4 * Sw_wSrw );
      aPs = aEs + aWs;
      sp (i) = aPs * S(i) - aEs * S(i+1) - aWs * S(i-1);

    }

    // Left boundary condition.

    /* */
    aP(bi) -= aW(bi);
    Sw_wSrw = S(bi-1) - Srw;
    Sw_wSrw4 = Sw_wSrw * Sw_wSrw * Sw_wSrw * Sw_wSrw;
    aWs = COEF4 * sqrt( Sw_wSrw4 * Sw_wSrw );  
    sp(bi) += ( 3.53e-6 - 3 * aWs * ( S(bi) - S(bi-1) ) );
    /* */

    // Update pressure at left boundary
    /* */
    phi(bi-1) = phi(bi)
      - aWs * (S(bi) - S(bi-1)) // / aW(bi) 
      + 3.53e-6 / (2); //* aW(bi));
    /* */

    aW(bi) = 0;

    // Right boundary condition.
    aP(ei) += aE(ei);
    Sw_eSrw = S(ei+1) - Srw;
    Sw_eSrw4 = Sw_eSrw * Sw_eSrw * Sw_eSrw * Sw_eSrw;
    aEs = COEF4 * sqrt( Sw_eSrw4 * Sw_eSrw );
    sp(ei) += ( 2 * aE(ei) * p(ei+1) +  aEs * ( S(ei) - S(ei+1) ) );

    aE(ei) = 0;

    return 0;
}

/*
 *  Lineal for realtive permeability , Upwind for Sw
 *
template<typename Tprec, int Dim>
inline bool WFIP1<Tprec, Dim>::calcCoefficients2D () 
{
  static prec_t Sw_e, Sw_w, S_n, S_s;
    static prec_t k_dx_mu_o = k / (dx * mu_o);
    static prec_t k_dx_mu_w = k / (dx * mu_w);
    static prec_t SrwSro1 = 1 - Srw - Sro;
    static prec_t SrwSro2 = SrwSro1 * SrwSro1;
    static prec_t SrwSro4 = SrwSro2 * SrwSro2;
    static prec_t SrwSro7 = SrwSro2 * SrwSro2 * SrwSro2 * SrwSro1;
    static prec_t COEF1 = k_dx_mu_w / SrwSro4;
    static prec_t COEF2 = k_dx_mu_o / SrwSro2;
    static prec_t COEF3 = k_dx_mu_o / ( SrwSro2 * SrwSro4 );
    static prec_t Sw_eSrw, Sw_eSrw4, Sw_wSrw, Sw_wSrw4;
    static prec_t SroSw_e2, SroSw_w2;
    static prec_t COEF4 = -k * p_threshold / 
      (mu_w * theta * dx * sqrt(SrwSro7));
    static prec_t aPs, aEs, aWs, aNs, aSs;

    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; sp = 0.0;


    for (int i =  bi; i <= ei; ++i) {  
      for (int j = bj; j <= ej; ++j) {  
	// Upwind
	if ( phi_0(i+1,j) >= phi_0(i,j) ) Sw_e = S(i+1,j);
	else                              Sw_e = S(i,j);
	if ( phi_0(i-1,j) >= phi_0(i,j) ) Sw_w = S(i-1,j);
	else                              Sw_w = S(i,j);

	if ( phi_0(i,j+1) >= phi_0(i,j) ) Sw_n = S(i, j+1);
	else                              Sw_n = S(i, j  );
	if ( phi_0(i,j-1) >= phi_0(i,j) ) Sw_s = S(i, j-1);
	else                              Sw_s = S(i, j  )
      
	Sw_eSrw  = Sw_e - Srw;
	Sw_eSrw4 = Sw_eSrw * Sw_eSrw * Sw_eSrw * Sw_eSrw;
	
	Sw_wSrw  = Sw_w - Srw;
	Sw_wSrw4 = Sw_wSrw * Sw_wSrw * Sw_wSrw * Sw_wSrw;
	
	SroSw_e2 = (1 - Sro - Sw_e) * (1 - Sro - Sw_e);
	SroSw_w2 = (1 - Sro - Sw_w) * (1 - Sro - Sw_w);
	
	aE (i,j) = COEF1 * Sw_eSrw4 + COEF2 * SroSw_e2 
	  - COEF3 * SroSw_e2 * Sw_eSrw4;
	aW (i,j) = COEF1 * Sw_wSrw4 + COEF2 * SroSw_w2 
	  - COEF3 * SroSw_w2 * Sw_wSrw4;
	aP (i,j) = aE (i,j) + aW (i,j) + aN (i,j) + aS(i,j);
      
	aEs = COEF4 * sqrt( Sw_eSrw4 * Sw_eSrw );
	aWs = COEF4 * sqrt( Sw_wSrw4 * Sw_wSrw );
	aPs = aEs + aWs;
	sp (i,j) = aPs * S(i,j) - aEs * S(i+1,j) - aWs * S(i-1,j);

    }

    // Left boundary condition.

    /* *
    aP(bi) -= aW(bi);
    Sw_wSrw = S(bi-1) - Srw;
    Sw_wSrw4 = Sw_wSrw * Sw_wSrw * Sw_wSrw * Sw_wSrw;
    aWs = COEF4 * sqrt( Sw_wSrw4 * Sw_wSrw );  
    sp(bi) += ( 3.53e-6 - 3 * aWs * ( S(bi) - S(bi-1) ) );
    /* */

    // Update pressure at left boundary
    /* *
    phi(bi-1) = phi(bi)
      - aWs * (S(bi) - S(bi-1)) // / aW(bi) 
      + 3.53e-6 / (2); //* aW(bi));
    /* *

    aW(bi) = 0;

    // Right boundary condition.
    aP(ei) += aE(ei);
    Sw_eSrw = S(ei+1) - Srw;
    Sw_eSrw4 = Sw_eSrw * Sw_eSrw * Sw_eSrw * Sw_eSrw;
    aEs = COEF4 * sqrt( Sw_eSrw4 * Sw_eSrw );
    sp(ei) += ( 2 * aE(ei) * p(ei+1) +  aEs * ( S(ei) - S(ei+1) ) );

    aE(ei) = 0;

    return 0;
}
    */

} // Tuna namespace


#endif //_WFIP1_H_

















